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We study instabilities and pattern formation in reaction-diffusion layers that are diffusively cou- 
pled. For two-layer systems of identical two-component reactions, we analyze the stability of ho- 
mogeneous steady states by exploiting the block symmetric structure of the linear problem. There 
are eight possible primary bifurcation scenarios, including a Turing- Turing bifurcation that involves 
two disparate length scales whose ratio may be tuned via the inter-layer coupling. For systems of 
n-component layers and non-identical layers, the linear problem's block form allows approximate 
decomposition into lower-dimensional linear problems if the coupling is sufficiently weak. As an 
example, we apply these results to a two-layer Brusselator system. The competing length scales en- 
gineered within the linear problem are readily apparent in numerical simulations of the full system. 
Selecting a \/2:l length scale ratio produces an unusual steady square pattern. 
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I. INTRODUCTION 

In 1952, Alan Turing hypothesized that reaction and 
diffusion could compete to create stationary spatial pat- 
terns [1]. This hypothetical mechanism for biological 
morphogenesis has been the theoretical foundation for 
decades of work on Turing patterns, which form when 
a rapidly diffusing activator interacts with a slowly dif- 
fusing inhibitor. Nearly 40 years later, experimentalists 
observed these patterns in a chemical reaction-diffusion 
system [2]. Since then, chemical systems have been the 
canonical testing ground for Turing patterns. 

A variation on the classic Turing system is the multi- 
layered system, in which each layer is a reaction-diffusion 
system that is diffusively coupled to adjacent layers. 
These coupled systems are common in the biological 
world, seen in neural, developmental, and ecological con- 
texts [3] . One example from neuroscience is a neural-glial 
network, consisting of a layer of neurons connected dif- 
fusively to a layer of glial cells, where each layer exhibits 
dynamics at different time scales. The chemicals released 
at a tripartite synapse (one glial cell and a pair of neu- 
rons) and their effect on those cells are known [4, 5]; 
however, the effect of glial cells on the network level re- 
mains a subject of ongoing study [6] . Understanding how 
coupled layers influence one another contributes to our 
understanding of these networks. 

Though experimental studies of the biological sys- 
tems are quite difficult, investigations of the fundamen- 
tal properties of coupled reaction-diffusion systems have 
progressed via chemical experiments. Experimentalists 
employ two thin gels (which contain the reactants) that 
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are put in contact with one another. By adding or re- 
moving a permeable membrane between the layers and 
by adjusting its properties, the coupling strength can be 
altered. This approach with the chlorine dioxide - iodine 
- malonic acid (CDIMA) reaction has produced superlat- 
tice patterns called black-eyed and white-eyed patterns 
which involve wavelength ratios of nearly 2:1; other ra- 
tios were not feasible for this reaction and experimental 
configuration [7] . Recent experiments have exploited the 
photosensitivity of the CDIMA chemical reaction, using 
an external light source to probe the interaction between 
different forced patterns [8]. For a broad overview of 
experimental and numerical results for some multi-layer 
systems, see [3]. 

A few theoretical studies of multilayer systems have 
taken place in the setting of diffusively coupled ordi- 
nary differential equations; this framework neglects spa- 
tial dependence within layers (and hence, spatiotcmpo- 
ral pattern formation) but is more easily analyzed than 
the spatial case. Linear stability analysis and numerical 
bifurcation studies reveal regimes of in-phase and out-of- 
phase oscillations of coupled Brusselators [9], as well as 
regimes of synchronization and chaos in coupled Orego- 
nators [10]. For Brusselators, regions of in-phase waves 
and echo waves, whose phase differs by half the period, 
can be determined analytically [11, 12]. 

Work incorporating spatial dependence within layers 
has also used linear stability and bifurcation analyses to 
determine and understand possible patterns, now in the 
setting of partial differential equation models of chemi- 
cal reactions. For coupled Oregonators, simulations re- 
veal twinkling eye patterns, Turing spots arranged in a 
hexagonal lattice that oscillate 120 degrees out of phase 
with their nearest neighbors, and traveling waves in Tur- 
ing structures, such as pinwheels in spots and traveling 
waves in labyrinths [13]. A numerically computed dis- 
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persion relation suggests that the twinkling eye pattern is 
due to an interaction of Turing and Hopf modes, whereas 
the traveling wave patterns are formed via a short wave 
instability [13]. Similar analyses have also elucidated the 
bifurcations to time-dependent Turing states and super- 
lattices in coupled Lengyel-Epstein equations as param- 
eters vary [14], in the presence of a delay [15], and with 
external forcing [16]. Simulations of coupled Brussela- 
tors demonstrate superposition patterns, twinkling eye 
patterns, and black-eyed and white-eyed superlattices. 
They are the result of two interacting Turing modes and 
occur when the ratio of the interacting modes is close to 
y/3:l, 2:1, 3:1, [17] or 4:1 [18]. The Jacobian matrix of 
this system has been studied numerically to understand 
these patterns [18, 19]. One analytical study of coupled 
Brusselators used a linear stability analysis to obtain con- 
ditions for existence of steady states and non-constant 
solutions [20]. Extending work on coupled layers, [21] 
studies networks of reaction-diffusion systems, which are 
closely related to the BZ-AOT experimental system [22]. 

Analytical calculations for layered reaction-diffusion 
systems can be difficult because of dimensionality. For 
instance, with m layers of n-componcnt reaction-diffusion 
systems, the linear problem is mnxmn; this suggests why 
even the linear results for the papers referenced above are 
largely numerical. In this paper, we show how the linear 
calculations may be simplified and harnessed to engineer 
certain aspects of nonlinear pattern formation. For the 
case of a two-layer, two-component system, we exploit 
the block symmetric form of the Jacobian to analyze the 
stability of homogeneous steady states. There are eight 
possible primary bifurcation scenarios, and we determine 
conditions under which each occurs. One possibility is 
a Turing- Turing bifurcation that involves two disparate 
length scales whose ratio may be tuned via the inter-layer 
coupling. For systems of n-component layers and non- 
identical layers, the linear problem's block form allows 
approximate decomposition into lower-dimensional linear 
problems if the coupling is sufficiently weak. We apply 
some results to a two-layer Brusselator system near the 
Turing- Turing bifurcation. The competing length scales 
engineered within the linear problem are readily apparent 
in numerical simulations of the full system. Selecting a 
\/2:l ratio produces a steady square pattern. Square su- 
perlattice Turing patterns have been previously reported, 
initially in [23], under the influence of external forcing. 
However, to our knowledge, a steady pattern of simple 
Turing squares (moreover, one obtained without forcing) 
has not been previously reported. 

The rest of this paper is organized as follows. Sec. II 
presents a linear stability analysis of coupled layers 
of reaction-diffusion systems, describing in detail the 
primary bifurcations for the case of two-layer, two- 
component systems. The linear algebra necessary to 
simplify the calculations is developed in the Appendices. 
Sec. Ill applies some of the results in order to engineer 
nonlinear patterns containing a desired length scale ratio, 
as demonstrated in simulations. We conclude in Sec. IV. 



II. LINEAR ANALYSIS OF COUPLED 
REACTION-DIFFUSION LAYERS 

We now present results for the (in)stability of trivial 
states of coupled reaction-diffusion layers. In Sec. II A 
through II C we focus on two-layer systems of identical 
two-component layers. Exploiting the block symmetric 
structure of the linearized problem, we find convenient 
expressions for the eigenvalues that are easily analyzed, 
and we enumerate the possible primary bifurcations. In 
Sec. II D, we mention a few brief results applying to 
systems with non- identical layers, more complicated cou- 
pling schemes, and systems with more chemical compo- 
nents. Because we begin with generic reaction-diffusion 
equations, the results are readily applied to specific sys- 
tems such as the Brusselator [24], the Lengyel-Epstein 
model [25], and so forth. 



A. Derivation of linearized problem 

We begin with nonlinear equations describing identical 
two-component reaction-diffusion layers that are coupled 
together, 

Ui =a(Uj - Ui) + F(Ui,Vi)+ V 2 U t , (la) 
V t = ${Vj - Vi) + G{U t , V l )+DV 2 V l . (lb) 

This model describes layers that are identical in their 
underlying chemical and physical properties. Through- 
out this section, i,j = 1,2, i ^ j indicates the layer. 
Ui(x, t), Vi(x, t) are chemical concentration fields, x is 
the spatial coordinate, t is time, and the over dot rep- 
resents a time derivative. The functions F and G are re- 
action kinetics terms whose functional form depends on 
the particular chemical model under consideration. The 
diffusivity of U is set to unity by a rescaling of the spa- 
tial coordinate; the diffusivity of V is D. Without loss 
of generality, assume V to be the more rapidly diffusing 
species, so that D > 1. Finally, a, f3 > are coefficients 
of diffusive coupling between the systems. 

We wish to study bifurcations from a spatially uni- 
form steady state. As pointed out in [14], different types 
of uniform states may be possible. One possibility is that 
the concentrations of the two layers are identical. A sec- 
ond possibility is that the two layers have distinct (uni- 
form) concentrations even though the underlying equa- 
tions are the same. In practice, the types of steady states 
that exist are determined by the particular form of the 
reaction kinetics functions F, G and the chemical param- 
eters therein. In this section, as the simplest case, assume 
that the two layers share the same uniform steady state. 
We relax this assumption in Sec. II D. 

Let the uniform steady state be Ui = U* , Vi = V* . 
Write the chemical fields as a perturbation around the 
steady state, and express the perturbation as a superpo- 



3 



sition of Fourier modes. 



Global extrema of trace and determinant 



u* 
v* 



E 



(2) 



The perturbation has wave number q = |q| and iij, q (i) 
and Vi^q(t) are Fourier wave amplitudes. The summation 
and the admissible q must be interpreted in a manner 
consistent with the boundary conditions of the govern- 
ing equations; for instance, if the equations are posed on 
an unbounded domain, then the summation is actually 
an integral over all q (per the Fourier transform). To as- 
sess the stability of the steady state, study the linearized 
problem governing the perturbations, 



q 2 Ui 



iii =a(uj — Ui) + aui+bvi 
ii = f3(vj — Vi) + cui +dvi — Dq 2 Vi 



(3a) 
(3b) 



For brevity, and as a convenient abuse of notation, we 
have suppressed the q dependence in the subscript of the 
Fourier wave amplitudes. The coefficients a, 6, c, and d 
are given by 



a = 



c = 



dF 
dU 
dG 
dU 



(U-,V) 



(U*,V*) 



b= 
d= 



dF 
dV 
dG 

W 



(U*,V) 



{U-,V) 



(4a) 
(4b) 



It is convenient to write the problem in matrix form. 
Let u = (ui, v\, u 2 , v 2 ) T . The linearized problem is 



Lu, 



P Q 
Q P 



L is of block symmetric form, with blocks 
P = 

Q = 



c 



a — q" — a b 

d-Dq 2 -(3 



a 

o p) ■ 



(5) 

(6a) 
(6b) 



We show in Appendix A that the eigenvalues of L are 
the eigenvalues of Li = P + Q and L 2 = P Q. Hence, 
the linear problem decomposes conveniently into two sub- 
problems described by the matrices 



a — q 

a - q 2 - 2a 



Jl ~ 1 c d-Dq 2 ^ 



b 

d - Dq 2 - 2/3 ,/ ' 



(7a) 
(7b) 



The matrix Li is simply the Jacobian corresponding to 
a solitary reaction diffusion layer. The effect of the cou- 
pling between layers is seen via L 2 . Though the full lin- 
ear problem is four-by-four with a quartic characteristic 
polynomial, the problem decomposes into these two 2x2 
problems, facilitating analysis. 



In Sec. IIC we will consider different bifurcation sce- 
narios by analyzing the trace Ti,2(<z) and determinant 
Ai i2 (g) of Li j2 , 



n(«) 


= a + d-(D + l)q 2 , 


(8a) 


T 2(<?) 


= T 1 -2(a + 0), 


(8b) 


Ai(<?) 


= Dq 4 - (aD + d)q 2 + ad - be, 


(8c) 


A 2 (g) 


= A 1 +2(aD + P)q 2 


(8d) 




+ 2(-ad- l3a + 2a/3). 





Here, we present two helpful observations. 

First, ti i2 (<7) are quadratic in q, each with a negative 
leading coefficient and no g 1 term, and hence have global 
maxima at q = 0. We have 



ri(O) = a + d, 

r 2 (0) = a + d-2(a + /3) <n{0). 



(9a) 
(9b) 



Second, Ai^(q) are even-powered quartics, each with 
a positive leading coefficient. Thus, these quantities 
have global minima. Label them (<?i, m i n , Ai iTO j„) and 
(g 2 jJn j n , A 2 tm i n ). Whether the global minima occur at 
zero or nonzero q depends on the sign of the quadratic 
coefficient. For Ai(q), 



UaD + d>0: 
2 aD + d 



Ai, r , 



2D ' 
(aD - d) 2 + ADbc 
4D ~ 



or 



If aD + d < 0: 



= o, 



= ad — be. 

Similarly, for A 2 (q), 

UaD + d- 2aD - 2/3 > 0: 
2 _aD + d- 2aD - 2/3 

^2, rain 



^2, rain 



2D 

(aD-d- 2aD + 2/3) 2 
\l) 



4Dbc 



or 



If aD + d - 2aD - 2/3 < 0: 



^2. rain 

A. 



= 0, 



^2,min = ad — be — 2ad — 2a/3 + 4a/3. 
Finally, note that q2, m in < Qi.min since a, /3 > 0. 



10a) 
10b) 



11a) 
lib) 



12a) 
12b) 



13a) 
13b) 
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TABLE I. Summary of possible primary bifurcations of the homogeneous steady state of (1). The four-dimensional linearized 
problem consists of two two-dimensional sub-problems per (7). We distinguish between two different classes of bifurcations. 
First, there are bifurcations due to eigenvalues in Li, which also occur in single-layer (traditional) two-component reaction- 
diffusion systems. These bifurcations are captured in Cases I - IV and are very well-known. Second, there are bifurcations due 
to eigenvalues in L 2 , and thus which depend on the diffusive coupling between the two layers. These are cases V - VIII. Below, 
a dash indicates no bifurcation, H indicates Hopf, T indicates Turing, and TH indicates Turing-Hopf. For each scenario, we 
state generic conditions on the traces and determinants Ti, 2 (g) and Ai, 2 (g) in (8). In practice, we enforce these conditions by 
controlling the global extrema of n,2(g) and Ai^(q); see Sec. II C for details. In the table, the wave number q c refers to a 
critical wave number; cases VII and VIII have two critical wave numbers. 

„ Bifurcation due to , , , . . , . . , . 

Case Ti(q) r 2 (g) Ai(g) A 2 (g) 



In L 



2 



I - - Ti(g)<0 r 2 (g)<0 Ai(qr)>0 A 2 (g) > 

II H n(0) = r 2 (g) < Ai(g) > A 2 (g) > 

n( 9 ^0) <0 

III T - Ti(q) < r 2 (g) < Ai(g c ) = A 2 (g) > 

Ai(g ^ q c ) > 

IV TH - n(0) = r 2 (g) < Ai(g c ) = A 2 (g) > 

n(g/0)<0 Ai(g/g c )>0 

V - T n(g)<0 r 2 (g)<0 Aj(g) > A 2 (g c ) = 

A 2 (g ^ q c ) > 

VI H T n(0) = r 2 (g)<0 Aj(g) > A 2 (g c ) = 

n(g/0)<0 A 2 (g/g c )>0 

VII T T n(0)<0 r 2 (g)<0 Ai( 9 i, c )=0 A 2 (g 2 , c ) = 

A 2 (g/gi, c )>0 A 2 (g/g 2 , c )>0 

VIII TH T n(0) = r 2 (g)<0 A^q^) = A 2 (g 2 , c ) = 

n(g/0)<0 A 2 (g/gi, c )>0 A 2 (g ± g 2 , c ) > 



C. Primary bifurcations 

We now consider possible primary bifurcation scenar- 
ios. Naively, Li and L 2 may each give rise to four dif- 
ferent primary bifurcation scenarios: none (linear stabil- 
ity), Hopf bifurcation (H), Turing bifurcation (T), and 
Turing-Hopf bifurcation (TH). Since the full linear prob- 
lem comprises Li i2 , there would be 4 x 4 = 16 primary 
bifurcation scenarios. 

However, due to the particular form of Li ;2 , any sce- 
nario involving a primary Hopf bifurcation in L 2 (that is, 
H or TH) is impossible. To see this, assume a primary 
Hopf bifurcation due to L 2 . This requires t 2 (0) = per 
(9). However, since t 2 (0) < ti(0) (with equality achieved 
only in the trivial case a = /? = 0) , the assumption means 
that a Hopf bifurcation would already have occurred due 
to Li, and hence the assumed bifurcation due to L 2 would 
not, in fact, be the primary one. Therefore, all bifurca- 



tion scenarios involving primary H or TH bifurcations 
due to L 2 are prohibited. This eliminates eight of the 
16 possible scenarios. Of course, if the layers were not 
identical, this result would not hold, and other primary 
bifurcations might be possible. For an example involving 
different Hopf bifurcations, see the nonspatial two-cell 
chemical model in [26]. 

The remaining eight possible primary bifurcation sce- 
narios are enumerated in Table I. We find the condi- 
tions for each case by analyzing Ti, 2 (q) and Ai ;2 (g) in 
the usual way to determine when a single eigenvalue or 
pair of eigenvalues crosses the imaginary axis, with all 
other eigenvalues contained in the left half of the com- 
plex plane. In these cases we distinguish between two 
different classes of bifurcations. First, there are bifurca- 
tions due to eigenvalues in Li, which also occur in single- 
layer (traditional) two-component reaction-diffusion sys- 
tems. These bifurcations are captured in Cases I - IV 
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and are very well-known. Second, there are bifurca- 
tions due to eigenvalues in L 2 , and thus which depend 
on the diffusive coupling between the two layers. These 
are Cases V - VIII. They correspond to Cases I - IV but 
with an additional Turing bifurcation due to L 2 . 

We apply the generic conditions in Table I to our spe- 
cific linear problem (5) by enforcing conditions on the 
global extrema of Ti^(q), Ai i2 (g). First, focus on the 
trace (the fourth and fifth columns of Table I.) An ex- 
amination of (8) and (9) shows that if Ti(0) < 0, then 
T i(q 7^ 0) < 0, and similarly for r 2 (q). Recall also, as 
noted in (9), that t 2 (0) < ti(0). Thus for our model, the 
condition that T%{q) < 0, required for all of the bifurca- 
tions in Table I, is subsumed in the condition on ti(0) 
and T\{q ^ 0) in that table. 

Now focus on conditions for the determinants (the 
sixth and seventh columns of Table I). These are eas- 
ily enforced by controlling Ai iTO j„ and A 2 , m i„ as given 
by (10) - (13). In cases of Turing bifurcations, the crit- 
ical WctvG numb cr s q \ c and/or q 2 . c arc identified with 
the locations of the global minima, namely q\ ;m i n and/or 

qi.min- 

There is still the matter of which expressions out of 
(10) - (13) apply for each case. For Cases I and II, 



cither (12) or (13) will apply for {q2, r , 



1), de- 



pending on chemical kinetics and parameters. If (12) 
applies, then (10) must apply for (gi >m i„, Ai :Jnin ) since 
qi,min < q\,min- If (13) applies, then one of (10) or (11) 
will apply for (q>i, TO ,„, Ai )min ), depending on chemical 
kinetics and parameters. Since a Turing bifurcation oc- 
curs at a nonzero wave number, (10) applies for Ai >m j n 
in Cases III and IV. In these cases, either (12) or (13) 
might apply for A 2iTOin , depending on chemical kinetics 
and parameters. Similarly, in Cases V - VIII, (12) ap- 



plies for A 2>r 

Ai rain- 



Since g 2>ri 



< 



qi,min, (10) applies for 



D. Extensions to other layered reaction-diffusion 

systems 

Suppose each layer comprises a reaction-diffusion sys- 
tem with n chemical components. Then generalizing (1), 
the governing equations are 



U, = Q(U, - Ui) + F(Ui) + DV 2 U,. 



(14) 



As before, i,j = 1,2. i ^ j indicates the layer. Uj(x, t) € 
R™ is a vector containing concentrations of the n chem- 
ical components in layer i. The vector function F £ R" 
describes reaction kinetics. The n x n diagonal matrix Q 
contains coupling coefficients, 



Q 







\ 







(15) 



a„J 



and the n x n diagonal matrix D contains diffusion coef- 
ficients, 

\ 



D 







V 







(16) 



Dn,J 



and V 2 is understood to operate on each element of Uj. 

Assume identical uniform steady states in each layer, 
Uj = U*. Then the linearized problem has the block 
structure (5), as in Sec. II A, only now 



P = dF 



u* 



q 2 Y> - Q, 



(17) 



and dF is the Jacobian of F. Of course, now P and 
Q are n x n matrices. Nonetheless, many features are 
preserved from the 2x2 case. The eigenvalues of the 
two-layer system still decompose into the eigenvalues of 



Li = P + Q = dF 
L 2 = P Q = dF 



u* 



-q 2 T>, 

- g 2 D - 2Q, 



(18) 
(19) 



where Li is simply the linear operator corresponding to 
a single (uncoupled) layer, and L 2 incorporates the effect 
of the coupling. 

Now, as in [14], allow the uniform steady state to com- 
prise different concentrations in each layer (even though 
the chemical parameters for each layer are identical) so 
that 



Ui=Ut, U 2 



Then the linearized problem is 



u = Lu, L = 




where 



dF 



7 2 D 



Q, 



S = dF 



u; 



q 2 D - Q. 



(20) 

(21) 

(22a) 
(22b) 



In this case, no simple formula exists for the eigenvalues 
of L in term of P, Q, and S. However, if we assume 
that coupling is weak, that is Q — > eQ where e <C 1 
then the eigenvalues of L are approximately equal to the 
eigenvalues of P and the eigenvalues of S. We show this 
in Appendix B. 

We may also suppose that the two layers have dis- 
tinct chemical kinetics. For instance, the system might 
be composed of two coupled Brusselators, but with a dif- 
ferent set of chemical control parameters selected for each 
layer. The governing equations for this case are 

Ui = Q(U 2 -U 1 ) + F 1 (U 1 )+D 1 V 2 U 1 , (23a) 

U 2 = Q(U 1 -U 2 )+F 2 (U 2 )+D 2 V 2 U 2 . (23b) 
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The two distinct chemical kinetics functions Fi ; 2 and the 
two distinct matrices of diffusion coefficients Di 2 reflect 
the different chemical parameters in each layer. The lin- 
earized problem has the same form (21), only now 



P = dFi 
S = dF 2 



u: 



9 2 Di Q, 
g 2 D 2 Q. 



(24a) 
(24b) 



The results of the previous paragraph still hold. For weak 
coupling, the eigenvalues are approximately those of P 
and S. 



III. MULTIPLE LENGTH SCALE SELECTION 

Sec. II showed that the uniform steady state of two 
identical, coupled reaction-diffusion layers may lose sta- 
bility via a codimension-two Turing- Turing bifurcation 
that involves two disparate wave numbers. We now ex- 
amine this bifurcation in more depth, and explore how 
the strength of coupling between the layers may be used 
to tune pattern selection and encourage the formation of 
spatial patterns with a desired length scale ratio. We 
apply results to the Brusselator in order to compute 
length scale ratios as a function of inter-layer coupling 
strength. Finally, we show via numerical simulation that 
we are able to engineer nonlinear patterns with pre- 
selected length scale ratios; this includes a steady square 
pattern. 



A. Length scale ratios 

We now focus on Case VII in Table I, which de- 
scribes the codimension-two Turing- Turing bifurcation. 
Our goal is to derive conditions for the Turing- Turing bi- 
furcation in terms of the parameters a, b, c, d, D, a, and 
(3, and to calculate the length scale ratio in terms of 
these parameters. Recall that the critical wave num- 
bers for a Turing- Turing bifurcation are qi_ c — q\, m in 
and Q2 lC = q_2,min as given by (10) and (12). 

The condition Ai jTO j„ = enforces a relationship be- 
tween a, b, c, d, and D, independent of the coupling 
parameters a and j3. The condition Ti(0) < means 
that a + d < 0. Therefore, a and d are oppositely signed. 
Recalling that (10) applies for Case VII, we know that 
aD + d > 0. In order for q\ c to be positive, a must be 
positive since D > 1. Hence, d < 0. For the remainder 
of this section, we assume that parameters satisfy these 
inequalities, 



a > 0, d<0, aD + d>0. 



(25) 



The next condition in Case VII is l^i,min = 0. Using 
(10b) and substituting (12b) yields 



from which two possibilities follow. Either 
(3 = aD - aD + d. 



or 



P = aD, 



(27) 



(28) 



The first possibility, (27), describes a line in a-(3 space, 
but the /3-intercept —aD + d is negative. Since a, /3 > 0, 
the condition is realizable only for 



a > 



aD -d 



D 



(29) 



Substituting (27), the wave number q 2 , c from (12) is 



<?2,c 



3aD -d- AaD 



2D 



(30) 



For a Turing bifurcation, q 2 _ c must be positive. Solving 
Q2.c > and (29) simultaneously leads to the inequal- 
ity a < 5d/D which cannot be satisfied because of (25). 
Hence, no Turing- Turing bifurcation is possible for (27). 

The second case, (28), also describes a line in a-(i 
space, but it emanates from the origin. Along this line, 
the wave number g 2 ,c is 



92, c 



aD + d- AaD 



2D 



which is positive so long as 
aD + d 



a < 



AD 



(31) 



(32) 



Thus, for < a < (aD + d)/AD and f3 = aD, 
codimension-two Turing- Turing bifurcations occur. The 
wave number ratio r q along this bifurcation curve is 



gl,c 



aD + d 



aD + d- 4aD ' 



(33) 



We will later use this result to choose chemical parame- 
ters giving rise to patterns dominated by a desired wave 
number (or alternatively, length scale) ratio. In an ex- 
periment, changing the coupling for two chemical species 
independently is generally not possible, and hence novel 
experimental approaches would be needed to fulfill con- 
dition (28). 

The issue of wave number ratios connects to resonant 
triad interactions, which are important to the study of 
some pattern selection problems. Our discussion here 
echoes in some respects the discussions of [27-30], which 
study resonant triads in Faraday waves. Resonant triad 
interactions, the lowest order nonlinear interactions, in- 
volve three modes with wave vectors Qi, Q 2 , and Q3 
satisfying the condition 



(aD - d) 2 = (aD-d- 2aD + 2/3) 2 , 



(26) 



Qi + Q 2 = Q 3 



(34) 
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FIG. 1. Diagram of resonant triads in Fourier space. The 
Fourier modes satisfy (34). Solid circles and vectors indicate 
neutral stability, and dotted ones indicate weak damping. In 
(a), | CJ 1,2 1 < IQ3 1 and the resonant angle satisfies < 9 res < 
2tt/3. In (b), |Qi, 2 | > |Q 3 | and 2tt/3 < 8 rea < ir. 



For the resonant triads that interest us, Qi ; 2 lie on a 
single critical circle in Fourier space, and Q3 is a weakly 
damped mode lying on a different, (nearly) critical circle. 
Eq. (34) determines an angle of resonance 9 res <G [0, 7r) 
between the two critical wave vectors via the trigonomet- 
ric relationship 



cos 



'res 

~2~ 



W3_ 

2Qi 



(35) 



If Qi < Q 3 
, s e (271-/3,71-). 



where |Qi| = |Q 2 | = Qi and |Q 3 | = Q 3 
then 9 res e [0,2tt/3). If Q x > Q 3 then 6, 
These two cases are pictured in Fig. 1. 

Resonant triad interactions may impact pattern se- 
lection. Heuristically, the interaction allows energy ex- 
change between the critical and damped modes. If the 
damped mode is a sink, drawing energy from the excited 
modes, the interaction is an anti-selection mechanism 
that suppresses patterns involving the resonant angle. 
Alternatively, if the damped mode is a source, feeding 
energy to the excited modes, patterns involving the res- 
onant angle - or equivalently, the associated length scale 
ratio - may be enhanced. 

For our reaction-diffusion system near the Turing- 
Turing bifurcation point, define Ai as the eigenvalue asso- 
ciated with qi iC having the largest real part; similarly for 
A 2 and g 2jC . Now detune slightly in parameter space from 
the Turing- Turing bifurcation, so that Ai j2 are small and 
oppositely signed. Consider the two different possibilities 
for resonant triads pictured in Fig. 1. First, assume that 
the critical modes have a smaller wave number than the 
weakly damped one, so that panel (a) applies. Recalling 
that q2. c < <Zi, c for the Turing- Turing bifurcation (we ex- 
clude the degenerate case of equality), this means that 
Qi = li.c an d Q 3 = q\ c . Combining (33) and (35) gives 
the resonance angle at the Turing- Turing point, 



cos 



'res 

IT 



1 aD + d 

2 V aD + d-AaD' 



(36) 



The right-hand side must be real and must not exceed 
unit magnitude. These requirements yield an admissible 



range of a in which our resonant triads exist, 
3 aD + d 



< a < 



16 D 



(37) 



which is a subset of the range in (32). 

For the alternate case in which the critical modes have 
a larger wave number than the weakly damped one, 
Fig. 1(b) applies. Then Qi = gi ;C and Q 3 = q 2tC , and 
the resonance angle is 



cos 



0, 



1 jaD + d-AaD 



2 J 2 V aD + d 
For this case, the entire range (32) is admissible. 



(38) 



B. Multiple length scales in coupled Brusselator 
layers 

As an example, we apply our results to the Brussela- 
tor [24]. For this chemical reaction, 



F(U,V) = A- (B + 1)U + U 2 V, 
G(U,V) = BU-U 2 V, 



(39a) 
(39b) 



in (1). A, B are chemical parameters. The steady state is 
(U*,V*) = (A, B/A) and the coefficients a,b,c,d in (7) 
are 

a = B-1, b = A 2 , c=-B, d=-A 2 . (40) 

For concreteness, take A = 3, B = 9, as do many of the 
examples in [17]. Then 



a = 8, 6 = 9, c=-9, d = -9. 



(41) 



To have Ai. min = in (10), the diffusion coefficient must 
be D = 2.25. Then 



<2i,c 



1.414. 



(42) 



To have a codimension two bifurcation that admits reso- 
nant triads, (28) must hold. Then 



<?2,c 



V2 - 2a. 



(43) 



For these chemical parameters, Fig. 2(a) shows the 
wave number ratio r q in (33) as a function of a at the 
Turing- Turing point. Fig. 2(b) shows the resonant triad 
angle 8 res in (36) and (38), also as a function of a. For 
the lower (solid) branch, the resonant triad corresponds 
to Fig. 1(a), in which the damped mode has larger wave 
number than the critical ones. For the upper (dashed) 
branch, the resonant triad corresponds to Fig. 1(b), in 
which the damped mode has smaller wave number. For a 
range of a, either branch is accessible, depending on how 
one detunes from the codimension-two point, i.e., which 
circle in Fourier space is damped. 
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FIG. 2. (a) Ratio r q in (33) of two (nearly) critical wave 
numbers near a Turing- Turing bifurcation in the Brusselator. 
The governing equations are (1) and (39) with A — 3, B — 9, 
D = 2.25. The coupling parameter a is a free parameter 
and P = aD. (b) Angle of triad resonance corresponding 
to (a). For the lower (dashed) branch, the resonant triad 
corresponds to Fig. 1(a), in which the damped mode has 
larger wave number than the critical ones. For the upper 
(dashed) branch, the resonant triad corresponds to Fig. 1(b), 
in which the damped mode has smaller wave number. See 
Sec. IIIB for details. 



C. Numerical simulation 

Using the linear stability results and the understanding 
of multiple critical length scales near the Turing- Turing 
bifurcation, we attempt to engineer patterns with desired 
ratios near the Turing- Turing point. As in Sec. IIIB, we 
adopt the Brusselator as our model and choose A = 3, 
B = 9 in (39). We pre-select a desired wave length ratio 
and set parameters to be very near the Turing- Turing 
bifurcation, but such that one of the (nearly) critical 
modes has maximum eigenvalue of 0.01 (and hence can 



grow) and the other (nearly) critical mode has maximum 
eigenvalue —0.01 (and hence is weakly damped). These 
conditions determine values of D, a, and j3. The compu- 
tational domain is periodic and square, with the length 
of each side eight times the wave length of the weakly 
growing mode. Starting from a random initial condition, 
we integrate the system in spectral space with 64 modes 
along each axis using the Expint exponential integrator 
package for Matlab [31] with a time step of h = 0.4 and 
Krogstad time-stepping. We run simulations to t = 4000, 
which for our parameter choices is long enough for the 
solution to approach an attractor. 

In our first example, we select the wave number ratio 
0.5sec(7r/12), corresponding to a resonant angle of 30°. 
These conditions determine D = 2.244, a = 0.723, and 
/3 = 1.633 (to three decimal places). Thus the damped 
mode has wave number q = 1.414, and the dominant 
mode has wave number q = 0.732. Fig. 3(a) visu- 
alizes this, showing the (analytically calculated) eigen- 
value with largest real part as a function of q. Fig. 3(b) 
shows the result of the full numerical simulation, namely 
a stripe-dominated pattern that is sometimes referred to 
as labyrinthine. The Fourier spectrum of this pattern 
in Fig. 3(d) shows active modes lying on two circles in 
Fourier space (though it is clearly not dominated by reso- 
nant triad interactions). From the radial power spectrum 
of the pattern in Fig. 3(c) (with units chosen so that the 
dominant peak is normalized to unity) we see that those 
circles correspond to the selected wave numbers. 

A more ambitious goal is to go beyond selecting a ratio 
of length scales and to actually select a particular pat- 
tern. In general, this requires nonlinear analysis. How- 
ever, we can show one example where harnessing the lin- 
ear stability results does lead to successful pattern selec- 
tion. For this case, we set the wave number ratio y/2:l so 
that the resonant angle is 90°. Optimistically, one might 
expect a square pattern, which is what we obtain in Fig. 
4(b). Fig. 4(d) shows that the angles between each dom- 
inant Fourier mode are 90 . The chemical parameters 
in this case are the same as the previous example, ex- 
cept that we have changed the coupling parameters to 
a = 0.490 and j3 = 1.113 in order to shift the (nearly) 
critical peak to the required value of gi jC w 1, shown in 
Figs. 4(a) and (c). Steady square patterns have been re- 
ported in photosensitive reaction diffusion systems forced 
with a square mask [32], and oscillatory square patterns 
have been observed in autonomous reaction-diffusion sys- 
tems with interacting Turing and Hopf modes [33]. We 
have not previously seen an unforced, steady square pat- 
tern reported in the chemical Turing pattern literature, 
and believe that our computational result in Fig. 4 rep- 
resents the first such example. 

To verify that the square pattern is robust to changes 
in domain size - and not dependent on having a com- 
putational domain whose side fits an integral number of 
wavelengths of the weakly growing mode - we repeat the 
calculation of Fig. 4 but use a box size of 5-\/3 w 8.7 
wavelengths per side rather than eight, as before. This 






■ 

* * 



FIG. 3. Numerical simulation of coupled Brusselators given by (1) with (39). Parameter values are A — 3, B — 9, D — 2.244, 
a — 0.723, and /3 = 1.633 (to three decimal places), (a) The eigenvalue with maximum real part is plotted as a function of 
wave number, (b) Striped pattern resulting from this choice of parameters. Dark and light regions indicated variations in 
concentration of chemical u in the top layer. The bottom layer looks the same but with light and dark regions reversed, (c) 
The radial power spectrum of the striped pattern with units chosen so that the dominant peak is normalized to unity, (d) The 
Fourier spectrum of the striped pattern. 



computation indeed still produces a square pattern, as 
shown in Fig. 5, albeit one with a different spatial orien- 
tation. 



IV. CONCLUSION 

Layered, spatially-extended reaction-diffusion systems 
are analytically taxing due to their (potentially high) di- 
mensionality. The intriguing laboratory experiments and 
numerical simulations of the past decade have been sup- 
ported by comparatively few theoretical works. In this 
paper, we have sought to develop some basic theory for 
simple layered scenarios, and to connect linear results to 
nonlinear pattern formation. 

First, we presented a linear stability analysis for cer- 
tain layered reaction-diffusion systems. For two-layer 
systems of identical two-component layers, we analyzed 
the stability of homogeneous steady states by exploiting 
the block symmetric structure of the linear problem. This 
analysis revealed eight possible primary bifurcation sce- 
narios, including a Turing- Turing bifurcation involving 
two length scales whose ratio may be tuned via the intcr- 



layer coupling. For systems of n-component layers and 
non-identical layers, the linear problem's block form al- 
lowed approximate decomposition into lower-dimensional 
linear problems for sufficiently weak coupling. 

We applied some results to a two-layer Brusselator sys- 
tem near the Turing- Turing bifurcation. We calculated 
the ratio of critical wave numbers as a function of the 
coupling parameter and harnessed the analytical results 
to pre-selected chemical and coupling parameters that 
should give rise to a particular ratio in a fully (weakly) 
nonlinear system. Numerical simulations indeed revealed 
patterns dominated by the chosen ratio. In one example, 
by pre-selecting a y/2:l ratio, we obtained (without exter- 
nal forcing of the system) a simple, steady square-lattice- 
based pattern. Our numerical simulations demonstrate 
potential applications of our results as a means of un- 
derstanding and engineering the instabilities in layered 
reaction-diffusion systems. However, to develop a more 
complete picture of pattern formation in these systems, 
nonlinear analysis is required. We expect future work 
could address these detailed questions of pattern selec- 
tion. 

Finally, we hope that our results might be of 
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Lb), 



-0.15 






(d) 


■ ■ ■ 




• • 




■ ■ ■ 





FIG. 4. Parameter values are A = 3, B — 9, D = 2.244, a — 0.490, and /3 = 1.113 (to three decimal places), (a) The eigenvalue 
with maximum real part is plotted as a function of wave number, (b) Square pattern resulting from this choice of parameters. 
Dark and light regions indicated variations in concentration of chemical u in the top layer. The bottom layer looks the same 
but with light and dark regions reversed, (c) The radial power spectrum of the square pattern with units chosen so that the 
dominant peak is normalized to unity, (d) The Fourier spectrum of the square pattern. 



use to experimentalists. For instance, the Lengyel- 
Epstein model of the two-layer CDIMA reaction [14] 
can be written with F(U, V) = A — U — AUV/{1 + U 2 ), 
G{U, V) = BU - BUV(1 + U 2 ) in (1). The coefficients 
a,b,c,d in (7) arc a = (3A 2 - 125)/7, b = -20A/j, 
c = 2A 2 B/j, and d = —5AB/j, where for convenience 
we define 7 = A 2 + 25. Assuming that the Turing- Turing 
bifurcation conditions of Sec. Ill A are met, the wave 
number ratio r q in (33) is \/2:l when 



3A 2 D - 5AB - 125D 
8D(A 2 + 25) 



(44) 



Experiments performed in this parameter regime might 
shed light on whether steady square-lattice-based pat- 
terns can indeed arise. 
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Appendix A: Eigenvalues of block matrices 

Here we show a useful identity for the eigenvalues of 
a block matrix with the symmetric form relevant to the 
stability calculation in Sec. II A - II C. 

First we perform a side calculation. Consider a block 
matrix of the form 



(Al) 




Assume S is invertible and factor this as 

1 q\ /p-qs-^ o\ 



L = 



S 



S !R 



(A2) 



where I is the (appropriately sized) identity matrix. Now 
apply results from [34] for determinants of block matrices. 
For the factors in (A2), we have 




dct(I)dct(S) =dct(S), 



(A3) 
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FIG. 5. Results analogous to Fig. 4, and with the same parameters. However, whereas the square computational domain in 
Fig. 4 had each side of length eight times the length of the weakly growing mode as determined from linear stability analysis, 
here we choose 5\/3 « 8.7 wavelengths per side in order to verify that the box size was not responsible for stabilizing the square 
pattern. Indeed, here we still obtain a square pattern, albeit one with a different orientation, (a) The eigenvalue with maximum 
real part is plotted as a function of wave number (identical to Fig. 4(a), reproduced here for convenience), (b) Square pattern. 
Dark and light regions indicated variations in concentration of chemical u in the top layer, (c) The radial power spectrum of 
the square pattern with units chosen so that the dominant peak is normalized to unity, (d) The Fourier spectrum of the square 
pattern. 



and 



sized square matrices. That is, 



det P-QS-R 
^ S *R I J 

= det(P - QS _1 R) det(I) 
= dct(P - QS X R). 

Combine (A2) - (A4) to obtain 

dct(L) = dct(S) det(P - QS _1 R). 



(A4a) 

(A4b) 
(A4c) 



(A5) 



We now turn to our main calculation of this Appendix. 
Consider the stability analysis in Sec. II A - II C, in which 
case R = Q, S = P in (Al) and P and Q are identically- 




(A6) 



Seek the eigenvalues by finding the roots of the character- 
istic polynomial Cl(A) = det(L — AI), or more explicitly, 



C L (A) =det 



P — AI Q 
O P -AI, 



(A7) 



Then Cl(A) takes the form 
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C L (A) = det[P-AI] det[P-AI-Q(P-AI) _1 Q], (A8a) 

= det[P - AI] 2 det[I - (P - AI) _1 Q(P - AI) _1 Q], (A8b) 

= det[P-AI] 2 dct[I-{(P-AI)- 1 Q} 2 ], (A8c) 

= dct[P-AI] 2 dct[I-(P-AI)- 1 Q]dct[I + (P-AI)- 1 Q], (A8d) 

= det[P- AI-Q]det[P- AI + Q], (A8e) 

= det[(P-Q)-AI]det[(P + Q)-AI], (A8f) 

= Cp_ Q (A)-C P+Q (A). (A8g) 



The first line follows from direct application of (A5). The 
second follows from pulling a factor of P — AI out of the 
second determinant and combining it with the first. The 
third line follows from noting the squared quantity. The 
fourth line follows from factoring a difference of squares. 
The fifth line follows from redistributing one factor of 
P — AI into each of the two other terms. The sixth 
line follows simply from commutativity of matrix addi- 
tion/subtraction, and the last line follows from the defi- 
nition of a characteristic polynomial. 



Appendix B: Eigenvalues of block matrices with 
blocks that are small in magnitude 

We now show an approximation for the eigenvalues for 
a block matrix of a particular form, where certain blocks 
are scaled by a small parameter. Begin with the matrix 



(Bl) 



which arises as the linearization of a problem considered 
in Sec. II D. In fact, P and S include additive factors of 
Q, so for convenience, we let P = P — Q and S = S — Q. 
For the case of weak chemical coupling, the entries in Q 
are small, so we let Q — > eQ where e < 1 is a small 
bookkeeping parameter. Our matrix now has the form 




Thus, the characteristic polynomial for (A6) factors 
into that of P — Q and P + Q, and therefore, the eigen- 
values of L in (A6) are the eigenvalues of P — Q and the 
eigenvalues of P + Q. 



L= /P-eQ eQ \ 
{ eQ S-eQj' 

The characteristic polynomial is 



(B2) 



C L (A) = det[S eQ AI] det[P eQ AI e 2 Q(S eQ AI) _1 Q] 

= det[S eQ AI]{det[P eQ AI] + C(e 2 )}, 

w det[S eQ AI] det[P eQ AI], 

- Cs(A)-C P (A). 



(B3) 
(B4) 

(B5) 
(B6) 



The first line follows from direct application of (A5). The 
second line follows from Jacobi's formula for the differ- 
ential of a determinant. The third line follows from ne- 
glecting the C(e 2 ) correction, and the final line follows 



from the definitions of S and P, and from the definition 
of a characteristic polynomial. Thus, the eigenvalues of 
(Bl) are approximately those of P and those of S so long 
as Q is scaled by a small parameter. 
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